//=========================================================================
//
// cutils.c
// 
// C-extensions for the implementation of the Gibson & Robinson model 
// in Python.
//
//=========================================================================

#include "Python.h"
#include "arrayobject.h"
#include <time.h>
#include "cutils.h"

#define firing(p, c)   ( p < c ? 1 : 0)
#define clip(value)    (value > 1 ? 1: value)

//=========================================================================
// 1. Set up the methods table                                           
//=========================================================================
static PyMethodDef _cutilsMethods[] = {
    {"randmatrix", randmatrix, METH_VARARGS},
    {"randvector", randvector, METH_VARARGS},
    {"clipped", clipped, METH_VARARGS},
    {NULL, NULL} /* sentinel, marks end of struc */
};

//=========================================================================
// 2. Initialize cutils functions                                        
//=========================================================================
void init_cutils() {
    (void) Py_InitModule("_cutils", _cutilsMethods);
    import_array(); // must be present for NumPY
}

//=========================================================================
/* 3. Functions code                                                     */
//=========================================================================
static PyObject *randmatrix(PyObject *self, PyObject *args){

    int  n;                    // size of the matrix 
    double c;                  // probability
    int seed;                  // seed
    npy_intp dims[2];
    PyArrayObject *matrix;      // output matrix
    int i,j;                    // to iterate rows and colums
    double myrand;              // random number

    // extract arguments
    if (!PyArg_ParseTuple(args, "idi", &n, &c, &seed)) return NULL;

    // set the seed to appropiate values
    if (seed == 0)
        srand(time(NULL));
    else
        srand(seed);

    // create NumPy array
    //dims[0] = dims[1] = n;
    dims[0] = dims[1] = n;
    //matrix = (PyArrayObject *)PyArray_FromDims(2, dims, NPY_INT);
    matrix = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_INT);

    // fill NumPy Array
    int * matrixdata = (int *)matrix->data;
    for(i=0; i<n; i++){
        for(j=0; j<n; j++){
            if (i==j)
                matrixdata[j + i*n] = 0;
            else{
                myrand = (double) rand() / RAND_MAX;
                matrixdata[j + i*n] = firing(myrand, c );
            }
        }
    }
    
    return PyArray_Return(matrix);
}

static PyObject *randvector(PyObject *self, PyObject *args){
    int  n;                     // vector dimension 
    double c;                   // probability
    int seed;                   // seed
    PyArrayObject *vector;      // array of vectors
    int dims[1];
    double myrand;              // random number
    int i;

    // extract arguments
    if (!PyArg_ParseTuple(args, "idi", &n, &c, &seed)) return NULL;

    // set the seed to appropiate values
    if (seed == 0)
        srand(time(NULL));
    else
        srand(seed);

    dims[0] = dims[1] = n;
    vector = (PyArrayObject *)PyArray_FromDims(1,dims, NPY_INT);
    //vector = (PyArrayObject *)PyArray_SimpleNew(1,dims, NPY_INT);

    // fill NumPy Array
    int * vectordata = (int *)vector->data;
    for(i=0; i<n; i++){
        myrand = (double) rand() / RAND_MAX;
        vectordata[i] = firing(myrand, c );
    }

    return PyArray_Return(vector);
}

static PyObject *clipped(PyObject *self, PyObject *args){

    PyArrayObject *Z;      // array of patterns 
    PyArrayObject *W;      // array of connections
    int * Zdata; 
    int * Wdata;
    int dims[1];
    int ncells, npatterns;
    PyArrayObject *output;      // array of vectors
    int * Jdata;
    int i, j;
    

    // extract arguments
    if (!PyArg_ParseTuple(args, "O!O!", 
        &PyArray_Type, &Z, &PyArray_Type, &W)) return NULL;

    Zdata = (int *)Z->data;
    Wdata = (int *)W->data;

    ncells = W->dimensions[1];
    npatterns = Z->dimensions[0];

    // create a matrix of zeros
    dims[0] = dims[1] = ncells; /* number of cells */
    //output = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_INT, 0);
    output = (PyArrayObject *)PyArray_FromDims(2, dims, NPY_INT);

    Jdata = (int *)output->data;

    // load synaptic weights
    // iterate through patterns
    for (i=0; i < npatterns; i++){ /* rows = number of patterns */

        for (j=0; j < ncells; j++){ /* cols = cell id */

            if (Zdata[j + i*ncells] == 1)
                Jdata[j + i*ncells] += Wdata[j + i*ncells]; 
        }
    }
    

    
    return PyArray_Return(output);
}
